#include "CGeometry.hpp"
#include "Option.hpp"

#include <cmath>
#include <string>
#include <sstream>
#include <fstream>
#include <stdexcept>
#include <algorithm>
#include <iterator>
#include <iomanip>
using namespace std;

CGeometry::CGeometry()
{
    point = NULL;
    elem = NULL;
}

CGeometry::~CGeometry()
{
    if(point != NULL)
    {
        for(size_t iPoint = 0; iPoint < _nPoint; iPoint++)
        {
            delete point[iPoint];
        }
        delete[] point;
    }

    if(elem != NULL)
    {
        for(size_t iElem= 0; iElem < _nElem; iElem++)
        {
            delete elem[iElem];
        }
        delete[] elem;
    }

    for(size_t iInterface = 0; iInterface< interface.size(); iInterface++)
    {
        if(interface[iInterface] != 0) delete interface[iInterface];
    }
}

void CGeometry::loadmesh(CConfig* config)
{
    string filename = config->getMeshFile();
    ifstream in;
    in.open(filename.c_str());
    if(in.is_open())
    {
        string textline;
        // load point
        while(getline(in, textline))
        {
            size_t position = textline.find("NDIME=",0);
            if (position != string::npos) 
            {
                textline.erase(0, 6);
                stringstream strline(textline);
                strline >> _nDim;
            }

            position = textline.find ("NELEM=",0);
            if (position != string::npos) 
            {
                textline.erase(0, 6);
                stringstream strline(textline);
                strline >> _nElem;
                elem = new CEntity*[_nElem];
            }

            position = textline.find ("NPOIN=",0);
            if (position != string::npos) 
            {
                textline.erase (0,6);
                stringstream strline(textline);
                strline >> _nPoint;
                point = new CEntity*[_nPoint];
                for(size_t iPoint = 0; iPoint < _nPoint; iPoint++)
                {
                    getline(in, textline);
                    stringstream pointline(textline);
                    double* coord = new double[_nDim];
                    size_t pointid = 0;
                    for(size_t iDim = 0; iDim < _nDim; iDim++)
                    {
                        pointline >> coord[iDim];
                    }
                    pointline >> pointid;
                    CEntity* tempPoint = new CEntity(pointid, _nDim);
                    tempPoint->setCG(coord);
                    point[iPoint] = tempPoint;
                }
                break;
            }

            position = textline.find ("NMARK=",0);
            if (position != string::npos) 
            {
                textline.erase(0, 6);
                stringstream strline(textline);
                long int nMarker;
                strline >> nMarker;
                config->setNMarker(nMarker);
            }
        }
        in.close();
         
        // load elem
        in.open(filename.c_str());
        while(getline(in, textline))
        {
            size_t position = textline.find ("NELEM=",0);
            if (position != string::npos) 
            {
                for(size_t iele = 0; iele < _nElem; iele++)
                {
                    getline(in, textline);
                    stringstream eleline(textline);
                    size_t eleType = 0;
                    eleline >> eleType;
                    size_t nTempElem = 0;
                    CEntity* tempElem;
                    if(eleType == TRIANGLE)
                    {
                        nTempElem = 3;
                    } else if (eleType == RECTANGLE)
                    {
                        nTempElem = 4;
                    }

                    size_t pointIndex;
                    CEntity** pointlist = new CEntity*[nTempElem];
                    for(size_t iPoint = 0; iPoint < nTempElem; iPoint++)
                    {
                        eleline >> pointIndex;
                        pointlist[iPoint] = point[pointIndex];
                    }

                    size_t eleId;
                    eleline >> eleId;



                    tempElem = NULL;
                    if(eleType == TRIANGLE)
                    {
                        tempElem = new CTriangle(eleId, _nDim);
                    } else if (eleType == RECTANGLE)
                    {
                        tempElem = new CRectangle(eleId, _nDim);
                    }
                    tempElem->init(pointlist);
                    elem[iele] = tempElem;
                }
                break;
            }

        }

        in.close();
        in.open(filename.c_str());
        // load bound
        while(getline(in, textline))
        {
            size_t position = textline.find ("NMARK=",0);
            if (position != string::npos) 
            {
                textline.erase(0, 6);
                stringstream strline(textline);
                long int nMarker;
                strline >> nMarker;
                config->setNMarker(nMarker);
                bound.resize(nMarker);

                for(long int iBound = 0; iBound < config->getNMarker(); iBound++)
                {
                    while(getline(in, textline))
                    {
                        size_t position = textline.find ("MARKER_TAG=",0);
                        if (position != string::npos) 
                        {
                            textline.erase(0, 11);
                            string markerTag;
                            stringstream strline(textline);
                            strline >> markerTag;
                            config->setMarkerTag(iBound, markerTag);

                            getline(in, textline);
                            position = textline.find ("MARKER_ELEMS=",0);
                            if (position != string::npos) 
                            {
                                textline.erase(0, 13);
                                long int markerNElem;
                                strline.clear();
                                strline.str(textline);
                                strline >> markerNElem;
                                config->setMarkerNElem(iBound, markerNElem);
                                for(long int iBoundElem = 0; iBoundElem < markerNElem; iBoundElem++)
                                {
                                    getline(in, textline);
                                    if(config->getNDim() == 2)
                                    {
                                        vector< long int > tempBound(3);
                                        strline.clear();
                                        strline.str(textline);
                                        strline >> tempBound[0] >> tempBound[1] >> tempBound[2];
                                        bound[iBound].push_back(tempBound);
                                    }
                                }
                            }
                            break;
                        }
                    }
                }
                break;
            }
        }




        in.close();

    } else
    {
        throw runtime_error("The in-file " + filename + " is not opened correctly!\n");
    }


}

void CGeometry::savemesh(CConfig* config)
{
    string filename = config->getOutMeshFile();
    ofstream out(filename.c_str());
    if(out.is_open())
    {
        out << "NDIME= " << _nDim << endl;
        out << "NELEM= " << _nElem << endl;
        for(size_t iele = 0; iele < _nElem; iele++)
        {
            out << setprecision(10) << elem[iele]->getType()<< "   ";
            for(size_t iPoint = 0; iPoint < elem[iele]->getNPoint(); iPoint++)
            {
                out << elem[iele]->getPoint(iPoint) <<  "   ";
            }
            out << elem[iele]->getId() << endl;
        }

        out << "NPOIN= " << _nPoint << endl;
        for(size_t iPoint = 0; iPoint < _nPoint; iPoint++)
        {
            for(size_t iDim = 0; iDim < _nDim; iDim++)
            {
                out << setiosflags(ios::scientific) << setprecision(10) << point[iPoint]->getCG(iDim) << "    ";
            }

            out << point[iPoint]->getId() << endl;

        }

        out << "NMARK= " << config->getNMarker() << endl;
        for(long int iMarker = 0; iMarker < config->getNMarker(); iMarker++)
        {
            out << "MARKER_TAG= " << config->getMarkerTag(iMarker) << endl;
            out << "MARKER_ELEMS= " << config->getMarkerNElem(iMarker) << endl;
            for(long int iBoundElem = 0; iBoundElem < config->getMarkerNElem(iMarker); iBoundElem++)
            {
                for(size_t iPoint = 0; iPoint < bound[iMarker][iBoundElem].size(); iPoint++)
                {
                    out << bound[iMarker][iBoundElem][iPoint] << "  ";
                }
                out << endl;
            }
        }
        out.close();
    } else
    {
        throw runtime_error("The out-file " + filename + " is not opened correctly!\n");
    }
}

void CGeometry::setIface(void)
{

    for(size_t iElem = 0; iElem < _nElem; iElem++)
    {
        for(size_t iface = 0; iface < elem[iElem]->getNFace(); iface++) 
        {
            size_t* pointlist = elem[iElem]->getFacePoint(iface);
            int aFace = findFace(pointlist, elem[iElem]->getNFacePoint());
            if(aFace == -1)
            {
                if(_nDim == 2)
                {
                    CEntity* tempIFace = new CLine(interface.size(), _nDim);
                    CEntity** aPointlist = new CEntity*[elem[iElem]->getNFacePoint()];
                    for(size_t iPoint = 0; iPoint < elem[iElem]->getNFacePoint(); iPoint++)
                    {
                        aPointlist[iPoint] = point[pointlist[iPoint]];
                    }

                    tempIFace->init(aPointlist);
                    tempIFace->parent.push_back(elem[iElem]->getId());
                    interface.push_back(tempIFace);

                    for(size_t iPoint = 0; iPoint < elem[iElem]->getNFacePoint(); iPoint++)
                    {
                        aPointlist[iPoint] = NULL;
                    }
                    delete[] aPointlist;
                    tempIFace = NULL;
                } else if(_nDim == 3)
                {
                }

            } else
            {
                interface[aFace]->parent.push_back(elem[iElem]->getId());
            }
        }
    }
}

int CGeometry::findFace(size_t* aPointList, size_t aNPoint)
{
    int tempFace = -1;
    for(size_t iface = 0; iface < interface.size(); iface++)
    {
        if(interface[iface]->getNPoint() == aNPoint)
        {
            size_t aPoint = 0;
            for(size_t iPoint = 0; iPoint < aNPoint; iPoint++)
            {
                if(point[interface[iface]->getPoint(iPoint)]->getId() == aPointList[0])
                {
                    aPoint = iPoint;
                    break;
                }
            }

            size_t tag = 0;
            for(size_t iPoint = 0; iPoint < aNPoint; iPoint++)
            {
                if (aPoint > interface[iface]->getNPoint())
                {
                    aPoint -= interface[iface]->getNPoint();
                }
                if(aPointList[iPoint] == interface[iface]->getPoint(aPoint))
                {
                    tag++;
                }
            }
            if(tag == interface[iface]->getNPoint())
            {
                tempFace = iface;
                break;
            }
        }
    }

    return tempFace;

}

void CGeometry::setPointNeighborPoint()
{
    for(size_t iElem = 0; iElem < _nElem; iElem++)
    {
        if(elem[iElem]->getType() == TRIANGLE || elem[iElem]->getType() == RECTANGLE)
        {
            for(size_t iPoint = 0; iPoint < elem[iElem]->getNPoint(); iPoint++)
            {
                //size_t* aPointlist = point[elem[iElem]->getPoint(iPoint)]->getFacePoint(iPoint);
                size_t* aPointlist = elem[iElem]->getFacePoint(iPoint);
                int flag = point[aPointlist[0]]->findNeighborPoint(aPointlist[1]);
                if(flag < 0)
                {
                    point[aPointlist[0]]->neighborPoint.push_back(aPointlist[1]);
                    point[aPointlist[1]]->neighborPoint.push_back(aPointlist[0]);

                    CEntity** tempPoint = new CEntity*[elem[iElem]->getNFacePoint()];
                    for(size_t jPoint = 0; jPoint < elem[iElem]->getNFacePoint(); jPoint++)
                    {
                        tempPoint[jPoint] = point[aPointlist[jPoint]];
                    }

                    CEntity* tempIFace = new CLine(interface.size(), _nDim);
                    tempIFace->init(tempPoint);
                    tempIFace->parent.push_back(elem[iElem]->getId());
                    interface.push_back(tempIFace);

                    point[aPointlist[0]]->neighborLine.push_back(tempIFace->getId());
                    point[aPointlist[1]]->neighborLine.push_back(tempIFace->getId());

                    delete[] tempPoint;
                    tempIFace = NULL;

                } else
                {
                    int lineTag = -1;
                    for(size_t iline = 0; iline < point[aPointlist[0]]->neighborLine.size(); iline++)
                    {
                        size_t templineTag = point[aPointlist[0]]->neighborLine[iline];
                        for(size_t jPoint = 0; jPoint < interface[templineTag]->getNPoint(); jPoint++)
                        {
                            if(aPointlist[1] == interface[templineTag]->getPoint(jPoint)) 
                            {
                                lineTag = templineTag;
                                break;
                            }
                        }

                    }

                    interface[lineTag]->parent.push_back(elem[iElem]->getId());
                }

            }
        }

    }
}

void CGeometry::setInterfaceNormal()
{
    for(size_t iFace = 0; iFace < getNIface(); iFace++)
    {
        if(interface[iFace]->getNParent() == 2)
        {
            size_t iPoint = interface[iFace]->getParent(0);
            size_t jPoint = interface[iFace]->getParent(1);
            if(_nDim == 2)
            {
                interface[iFace]->normal = new double[2];
                double tempx = interface[iFace]->getCG(0);
                double tempy = interface[iFace]->getCG(1);
                double coord1[2];
                coord1[0] = point[interface[iFace]->getPoint(0)]->getCG(0);
                coord1[1] = point[interface[iFace]->getPoint(0)]->getCG(1);

                double coord2[2];
                coord2[0] = point[interface[iFace]->getPoint(1)]->getCG(0);
                coord2[1] = point[interface[iFace]->getPoint(1)]->getCG(1);
                double tempNormal[2];
                tempNormal[0] = coord2[1] - coord1[1];
                tempNormal[1] = coord1[0] - coord2[0];

                double temp = tempNormal[0] * (tempx - elem[iPoint]->getCG(0)) 
                    + tempNormal[1] * (tempy - elem[iPoint]->getCG(1));
                if (temp < 0.0)
                {
                    tempNormal[0] *= -1;
                    tempNormal[1] *= -1;
                }
                if(iPoint < jPoint)
                {
                    interface[iFace]->normal[0] = tempNormal[0];
                    interface[iFace]->normal[1] = tempNormal[1];
                } else
                {
                    interface[iFace]->normal[0] = -tempNormal[0];
                    interface[iFace]->normal[1] = -tempNormal[1];
                }
            }
        } else if(interface[iFace]->getNParent() == 1)
        {
            size_t iPoint = interface[iFace]->getParent(0);
            if(_nDim == 2)
            {
                interface[iFace]->normal = new double[2];
                double tempx = interface[iFace]->getCG(0);
                double tempy = interface[iFace]->getCG(1);
                double coord1[2];
                coord1[0] = point[interface[iFace]->getPoint(0)]->getCG(0);
                coord1[1] = point[interface[iFace]->getPoint(0)]->getCG(1);

                double coord2[2];
                coord2[0] = point[interface[iFace]->getPoint(1)]->getCG(0);
                coord2[1] = point[interface[iFace]->getPoint(1)]->getCG(1);
                double tempNormal[2];
                tempNormal[0] = coord2[1] - coord1[1];
                tempNormal[1] = coord1[0] - coord2[0];

                double temp = tempNormal[0] * (tempx - elem[iPoint]->getCG(0)) 
                    + tempNormal[1] * (tempy - elem[iPoint]->getCG(1));
                if (temp < 0.0)
                {
                    tempNormal[0] *= -1;
                    tempNormal[1] *= -1;
                }
                interface[iFace]->normal[0] = tempNormal[0];
                interface[iFace]->normal[1] = tempNormal[1];
            }
        }
    }
}

void CGeometry::setElemVolume()
{
    for(size_t iElem = 0; iElem < _nElem; iElem++)
    {
        size_t type = elem[iElem]->getType();
        double tempVolume = 0.0;
        switch (type)
        {
            case TRIANGLE:
                double coordTri[3][2];
                for(size_t iPoint = 0; iPoint < 3; iPoint++)
                {
                    coordTri[iPoint][0] = point[elem[iElem]->getPoint(iPoint)]->getCG(0);
                    coordTri[iPoint][1] = point[elem[iElem]->getPoint(iPoint)]->getCG(1);
                }
                tempVolume = 0.5 * (
                        (coordTri[0][0] - coordTri[1][0]) * (coordTri[2][1] - coordTri[1][1])
                        - (coordTri[2][0] - coordTri[1][0]) * (coordTri[0][1] - coordTri[1][1])
                        );

                elem[iElem]->setVolume(fabs(tempVolume));
                break;
            case RECTANGLE:
                double coordRec[4][2];
                for(size_t iPoint = 0; iPoint < 4; iPoint++)
                {
                    coordRec[iPoint][0] = point[elem[iElem]->getPoint(iPoint)]->getCG(0);
                    coordRec[iPoint][1] = point[elem[iElem]->getPoint(iPoint)]->getCG(1);
                }
                tempVolume = 0.5 * ( (coordRec[0][0] - coordRec[2][0]) * (coordRec[1][1] - coordRec[3][1])
                        + (coordRec[3][0] - coordRec[1][0]) * (coordRec[0][1] - coordRec[2][1]));

                elem[iElem]->setVolume(fabs(tempVolume));
                break;
        }
    }
}

void CGeometry::setBound(CConfig* config)
{
    long int nMarker = config->getNMarker();
    for(long int iMarker = 0; iMarker < nMarker; iMarker++)
    {
        long int nBound = config->getMarkerNElem(iMarker);
        for(long int iBound = 0; iBound < nBound; iBound++)
        {
            CEntity** tempPointList = new CEntity*[2];
            tempPointList[0] = new CPoint(bound[iMarker][iBound][1], 2);
            tempPointList[1] = new CPoint(bound[iMarker][iBound][2], 2);

            CEntity* templine = new CLine(0, 2);
            templine->init(tempPointList);

            bool isBound = false;
            for(size_t iInterface = 0; iInterface < interface.size(); iInterface++)
            {
                isBound = checkSameElem(templine, interface[iInterface]);
                if(isBound) 
                {
                    long int tempBoundType = config->getMarkerType(iMarker);
                    interface[iInterface]->setBoundType(tempBoundType);
                    break;
                }
            }
            
            
            delete tempPointList[0];
            delete tempPointList[1];
            delete[] tempPointList;

            delete templine;

        }

    }
}

bool CGeometry::checkSameElem(CEntity* aElem, CEntity* bElem)
{
    bool flag = false;
    long int aPoint0 = aElem->getPoint(0);
    long int aPoint1 = aElem->getPoint(1);

    long int bPoint0 = bElem->getPoint(0);
    long int bPoint1 = bElem->getPoint(1);

    if((aPoint0 == bPoint0 && aPoint1 == bPoint1) 
            || (aPoint0 == bPoint1 && aPoint1 == bPoint0))
        flag = true;
    return flag; 
}
